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solve the time-dependent Schrodinger equation for many-electron atoms, and particularly apply it 
to the treatment of photoionization processes in atoms. The method is presented in a very general 
formulation and incorporates a wide range of commonly used approximation schemes, like the single- 
active electron approximation, time-dependent Configuration Interaction with single-excitations, or 
the time-dependent R-matrix method. We proof the applicability of the method by calculating the 
photoionization cross sections of Helium and Beryllium. 
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I. INTRODUCTION 

The numerical simulation of quantum mechanical 
many-body systems is hampered by the exponentially 
growing effort required to directly solve the time- 
dependent Schrodinger equation (TDSE). As a conse- 
quence, only one-particle systems can be described in 
full generality. For two-particle systems (and reduced 
two-particle systems) — although great progress has been 
made in the last decade, e.g. in the description of dou- 
ble ionization of Helium [U |2] , molecular Hydrogen [3 or 
Beryllium [IHn! — there exist several scenarios for which 
direct solutions are still hardly or even not feasible: 
among them is, e.g., pump/probe spectroscopy [311], the 
laser assisted Auger decay [9j [10] , and correlated high- 
order harmonic generation [111 112) , all of which require 
large angular momentum expansions to obtain converged 
results. Further, the first-principle treatment of few- and 
many-electron systems is becoming more and more im- 
portant. This can be seen, for instance, in the measure- 
ments of a time-delay in the photoionization from differ- 
ent atomic shells of Neon and Argon atoms [13] , which up 
to now could not been reproduced in simulations. Other 
important processes which would profit from direct so- 
lutions, beside those already mentioned, are tunnel ion- 
ization [M] and excitation of hole states [HI [16] . In this 
work, we introduce a method to the field of numerical 
simulations of photoionization, which is capable of treat- 
ing a broad class of time-dependent physical processes, 
including the ones mentioned, in a very general manner. 

Several methods have been designed to attack the 
many-particle TDSE by using a reduced description, 
such as, e.g., time-dependent density functional the- 
ory [17], Nonequilibrium Green functions [I8H20] or 
semi-empirical approaches [ST] [55]; they often are ap- 
plicable to rather large systems, but typically lack a 
control of the accuracy. Regarding the wavefunction 
based schemes, we mention the Multiconfiguration time- 
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dependent Hartree-Fock (MCTDHF) method [Ml [H- 
[26], the time-dependent Configuration Interaction sin- 
gles method (TD-CIS) P7| and the single- and two- 
active electron approximation (SAE/TAE) [Mllin]- Fur- 
ther, there is the time-dependent R-matrix method 
(TDRM) [30l [31] , which can be considered the most suc- 
cessful time-dependent approach to few-electron systems 
so far. The restricted active space method employed in 
this work can be considered as a superset to all the men- 
tioned determinant schemes: it contains the TD-CIS, 
SAE/TAE methods, and even TDRM as special cases, 
and it can also be used to extend the applicability of the 
MCTDHF method. As will be discussed later, due to its 
generality, we believe that it should be applicable to pro- 
cesses which up to now could not be efficiently described 
in direct calculations. 

In this work, we aim at presenting the first test cal- 
culations, without fully exploiting the capabilities of the 
TD-RASCI method. Rather, as a first step, we focus 
on a standard problem in photoionization, namely the 
calculation of total cross sections of the atoms Helium 
and Beryllium. We therefore employ an explicitly time- 
dependent description, and compare different approxima- 
tions for the wavefunction with experimental and theo- 
retical results from the literature. Our results show that 
the TD-RASCI method yields an accurate description of 
the doubly-excited states at a fraction of the effort of full 
direct solutions. 

The paper is organized as follows: section II introduces 
the full and restricted configuration interaction schemes, 
which are subsequently related to commonly used ap- 
proximations in the treatment of photoionization, and 
gives the main ideas of our numerical implementation. In 
section III, we particularly focus on the treatment of pho- 
toionization processes. Therefore, we introduce a parti- 
tioning of the coordinate space, as well as a mixed single- 
particle basis, which will turn out a very convenient in- 
gredient of the present method. Section IV presents the 
first numerical results for Helium and Beryllium. In sec- 
tion V, we summarize the concepts and results. 
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II. TIME-DEPENDENT RESTRICTED ACTIVE 
SPACE CONFIGURATION INTERACTION 
(TD-RASCI) 

A. Full Configuration Interaction 

Our aim is to solve the time-dependent Schrodinger 
equation (TDSE) 

tdt\^{t)) = H{t)\^{t)), (1) 

for the iV-particle wavefunction | vE'(i)), which provides 
the complete information of the system (in a pure state 
description). We concentrate on a Hamiltonian describ- 
ing Coulomb-interacting fermions in an atom of charge Z 
subjected to an external laser field (we use atomic units) , 

(2) 

The method presented in the following, however, is com- 
pletely general and might be applied to a large variety of 
other physical systems. We begin with the expansion of 
the wavefunction in a set of Slater determinants, 

= E^HO (3) 

I en 

in which the multi-index / = (ii, • • • , i^) e Q specifies 
the occupied single-particle spin-orbitals | V'fc ) i which are 
assumed to be orthonormal throughout. The sum is per- 
formed over an index-set ft C N^, which determines the 
set of Slater determinants included in the expansion, and 
thus the accessible subspace of the Hilbert space Hn in 
which the wavefunction lives. By insertion of the ansatz 
([3]) into the TDSE, we obtain the equation of motion for 
the expansion coefficients, 

iCi{t) = ^(/|i7(i)| j)C,7(t), (4) 
Jen 

which is just the TDSE projected onto the subspace of 
T-Ln defined by fl. In order to solve this equation, the 
arising Hamiltonian matrix elements have to be evaluated 
using Slater-Condon rules, see e.g. Ref. [32]. 

To proceed, one needs to specify the index set 57. By 
choosing 

f^FCI = {(«!,• ■• ,iAr) I 1 < «1 < ■■ • < «Ar < '^N^,] , (5) 

one obtains the Full Configuration Interaction (FCI) 
ansatz. The determinant basis has a size of (2^") and 
spans the maximal accessible Hilbert space for a given 
spin-orbital basis of size 2Nb. The corresponding wave- 
function and operators are, up to the discretization and 
errors in the time-integration, represented exactly. In 
particular, the wavefunction contains the entire possi- 
ble correlation, and its solution provides the benchmark 



result for any approximate method in the same single- 
particle basis. In the limit of a sufficiently accurate 
single-particle basis, one essentially recovers the true re- 
sult. 

Eqs. ^ and Q present the basis of all direct ap- 
proaches, though they are often stated in a different form. 
In the photoionization community, and particularly in 
the treatment of two-particle systems, one often encoun- 
ters the close-coupling ansatz, where the expansion is 
made in angular-momentum eigenfunctions. In quantum 
chemistry, on the other hand, it is common to work in a 
basis of spin-eigenfunctions, so called configuration state 
functions [SS] . Both types of basis sets may be obtained 
from the Slater-determinant basis by taking appropriate 
linear combinations, the expansion coefficients of which 
are given in terms of Clebsch-Gordan coefficients. In a 
symmetry adapted basis set, one typically has to deal 
with less basis states and also obtains a clearer inter- 
pretation of the wavefunction, but this often goes at the 
cost of a more difficult evaluation of Hamiltonian matrix 
elements (particularly for more than two particles). 

As mentioned previously, the time-dependent Full 
Configuration Interaction method can hardly be applied 
to photoionization processes of TV > 2 particle systems 
(in three dimensions) , since the modelling of an adequate 
continuum often requires a rather large single-particle ba- 
sis. This causes the FCI expansion to become unfeasibly 
large. A trivial solution is to employ smaller basis sets 
resp. grids, as it is done in quantum chemistry calcula- 
tions or in a recent time-dependent close-coupling study 
of the Lithium atom |33]. The RAS method presented 
in the next section takes another approach: it retains 
the accuracy of the single-particle basis, but restricts the 
wavefunction on the many-body level. 



B. Restricted active space Configuration 
Interaction 

Time-dependent restricted active space Configuration 
Interaction (TD-RASCI) provides a way to effectively re- 
duce the determinant basis size. The underlying idea is 
quite simple: given the huge Full-CI space, remove all 
parts of this space which expectedly will not be occu- 
pied by the wave function. Practically, this means that 
a certain set of Slater determinants is dropped from the 
expansion (|3|, which reduces the size of the discretized 
Hilbert space and thus facilitates the numerical solution. 
The RAS method formalizes this idea and provides a sys- 
tematic approach for the selection of the important de- 
terminants. It is routinely used in quantum chemistry for 
more than 20 years ^35^, and there especially in Config- 
uration Interaction and Multiconfiguration Hartree-Fock 
calculations. In this context, it is also termed restricted 
active space self-consistent field (RASSCF). The idea 
of selecting individual configurations, however, has al- 
ready been used much earlier and was applied, e.g., in 
Configuration Interaction singles- and/or doubles calcu- 
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FIG. 1. Illustration of the restricted active space scheme 
for a number of P = 4 partitions of the single-particle basis 
B. In each partition, one imposes certain restrictions on the 
allowed particle numbers. The total A^-particle wavefunction 
I 4' ) is constructed as the tensor product of the A^i-particle 
wavefunctions j ), with Ni — N. 



lations In the time-dependent treatment of pho- 

toionization processes, up to now, we are only aware of 
special cases, such as the single-active electron approxi- 
mation (SAE) or time-dependent Configuration Interac- 
tion singles (TD-CIS). They will be related to the present 
method later in this work. 

The RAS method is formally obtained by imposing 
restrictions on the set Q of allowed Slater determinant 
indices. This, however, can be cumbersome without a 
clear and intuitive picture in mind. Therefore, a descrip- 
tive procedure may be obtained with the following steps: 



(i) Divide the single-particle basis 

13 {| -01 ),•••, I ■02Wi, )} 
into an arbitrary number P of partitions Bi, 



with 



B ^ 6i U • • • U 6^ 



(6) 



(7) 



(8) 



A division is thus defined by the P + 1 numbers 
(Po = 1, Pi, • • ■ ,PP-i, PP = 2iVb), i.e. by P - 1 
free parameters. The number of orbitals in the par- 
tition Bi is denoted by N^^i] by definition, one has 
J2j ^b,j ~ 2Nij. The partitions are visualized by 
the four black boxes in Fig. [T] 

(ii) Impose restrictions on the allowed particle num- 
bers. Therefore, for each partition Bi, we specify 
the minimal and maximal particle number, A^min,i 
and A'max i , and allow only for particle numbers Ni 



The restrictions should be assigned to match the oc- 
curring physical processes as good as possible, but 
at the same time result in only a moderate number 
of determinants. 

Each pair (Bi, Ni) obtained this way is related to a dis- 
crete Hilbert space Hi{Ni), which is the span of all A^^- 
particle Slater determinants constructed from the trun- 
cated single-particle basis Bi. For an example, see the 
right-hand part of Fig. [T] The total Hilbert space Hras 
is thus decomposed as 

Hras - U Hi(Ari) X ••• xHp(iVp), (9) 

i.e. as the unification of the cartesian products of all 
sub-Hilbert spaces which have the correct particle num- 
ber N and satisfy the RAS constraints. It has the total 
dimension 



dim('H 



RAS ) — 



E 



Ni 



Nb,P 
Np 



(10) 



JV„i„.i<Afi<Ar,„ 



and the wavefunction is set up as the antisymmetrized 
product of the given Slater determinants (which is just a 
determinant in the original Hilbert space), 



,{t) A\h)---\lp). (11) 



Finally, in order to bring the wavefunction expansion to 
the form of Eq. ([s]), we need to specify the set of allowed 
Slater determinant indices, 



RAS 



{{ii,- ■■ ,iN)€ ^FCi I Vje{l,--- ,P} : (12) 



Nnunj < Nj < Nn,in,j } , 



Here, the occupation Nj of the j—th partition Bj is given 
by 



N . 

k=l 



1, 

0, 



Pj-1 < ik < Pj 
else 



(13) 



Note that the notation ft^iAS does not explicitly state 
the strong dependence on the RAS parameters, i.e. on 
the partitioning {pi, ■ ■ ■ ,pp_i) and the range of allowed 
particle numbers defined by A^min.i and A'max.i- 

If we would have made no restriction on the particle 
numbers, i.e. allowed for < A^i < A^ in each partition 
Bi, we would essentially recover the Full Configuration 
Interaction method. By restricting the accessible many- 
body Hilbert space in the way just presented, one can 
significantly reduce its dimension and thus enable time- 
dependent Configuration Interaction calculations which 
are far beyond reach of the Full CI scheme. 
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FIG. 2. Special cases of the TD-RASCI scheme for the exam- 
ple of Beryllium. The numbers Ni label the allowed particle 
numbers in the partitions Bi (black boxes). Active electrons 
are shown as open (red) circles. Left: Single-active electron 
approximation with an active 2s orbital. The gray-shaded or- 
bitals are fixed, and no more than a single electron is allowed 
in the continuum. Right: Time-dependent Configuration In- 
teractions singles. All electrons are active, but only single 
excitations from the groundstate are included. 



C. RAS examples: SAE, TAB, TD-CIS 

The previous section might be regarded a bit technical, 
as we considered the formalism in a very general fash- 
ion. Actual applications of the RAS scheme are often 
much simpler. To get a feeling for its capabilities, the 
present section shall give the first examples where the 
RAS scheme reduces to (variants of) some well-known 
methods in the treatment of photoionization. At the 
same time, we wish to show how natural one can ar- 
rive at an extension of these standard methods. Let 
us start with the single-active electron (SAE) approxi- 
mation. There, the A'^-particle problem is reduced to a 
single-particle problem by freezing all electrons except 
the active one. An example of a corresponding basis 
partitioning for Beryllium is shown in the left part of 
Fig. [2j The total number of Slater determinants is given 
by 2Ni, — N, i.e. one determinant for each spin-orbital 
the active electron can reach. If we further include the 
spin-symmetry requirement (Sz) — 0, the number of de- 
terminants reduces to A^;, — N, as then only the spin-up 
orbitals are accessible. Similarly, also (Lz) = could be 
imposed. With the effort scaling only linearly with Ni,, 
the solution of the corresponding equation of motion ^ 
has only the complexity of a single-particle problem. 

In the particular form commonly applied in strong-field 
physics, the SAE approximation has been criticized for 
not being an ab-initio method |37j . This is due to the re- 
duction of an A^-particle problem to merely a one-particle 
equation, and the related need to estimate the effective 
potential experienced by the active electron. We want to 
stress, that this does not hold for the SAE approximation 
presented here, as the whole derivation assumed through- 
out an A^-particle representation. Similar variants of this 
extended SAE approximation have already been given in 



Refs. [371 13H]- Moreover, the scheme is easily extendable: 
in the case that another electron is likely to get ionized, 
one could treat the related orbital active as well, and 
thereby obtain the two-active electron (TAE) approxi- 
mation. The TAE approximation has been applied, for 
instance, to Beryllium [S] |6l [29|. In these works, how- 
ever, like in the SAE approximation one has to assume 
effective potentials, which are often of Hartree-Fock type 
and held fixed during the propagation. In contrast, the 
improved TAE approximation within the RASCI scheme 
contains the exact interaction terms, which are consis- 
tently adjusted during the propagation. The effective 
two-particle problem obtained in the TAE approxima- 
tion can already be difficult to solve for large basis sizes 
Ni,. Therefore, in order to retain the single-particle com- 
plexity, one could further restrict the number of allowed 
electrons in the virtual orbital space to one, which would 
lead to an only twice as large problem compared to a 
single active electron. 

Repeated application of this idea leads to the time- 
dependent Configuration Interaction singles (TD-CIS) 
approximation, where each of the N electrons is consid- 
ered as active, but only a single particle is allowed in 
the virtual orbital space at the same time. The visual- 
ization of this constraint is shown in Fig[2j The TD-CIS 
scheme has been applied with success to processes such as 
high-order harmonic generation in noble gas atoms [IT] 
or hole excitation processes [15]. It naturally lacks, how- 
ever, the description of transitions leading to doubly and 
higher excited states. As will be shown later in the re- 
sults section, the RAS scheme makes it easy to add cer- 
tain selected states, and by this can account e.g. for 
the doubly-excited resonances in Helium. The TD-CIS 
approximation could also be extended to include all dou- 
ble (TD-CISD) and higher (TD-CISDT, etc) excitations, 
however, at the price of a significantly increased effort. 



D. Numerical implementation 

In the following we summarize the main ideas of our 
numerical implementation of the TD-RASCI method. 
One of its great advantages is that, when properly im- 
plemented, it comprises a lot of different approximations 
in a single and generic program, which are achieved by a 
simple change of the RAS parameters. Restricted active 
space CI follows the basic work cycle as found in most CI 
implementations. In particular, the steps to accomplish 
are 

(1) Choose the number of particles N, and an appropri- 
ate single-particle basis B = {Itpk}}- 

(2) Select the RAS constraints, i.e. a partitioning of the 
basis into P elements Bi and corresponding restric- 
tions on the particle numbers Ni. 

(3) Construct the determinant basis and calculate the 
Hamiltonian matrix using the Slater-Condon rules. 
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(4) Solve the time-independent Schrodinger equation to 
obtain the initial eigenstate. 

(5) Integrate the time-dependent Schrodinger equation 
up to the required system time T. 

The first step is similar to Full CI calculations, al- 
though some care has to be devoted to the selection of the 
single-particle basis. Our choice is described in detail be- 
low, in section |III C[ essentially, the basis should be flex- 
ible enough to describe each possible state the electrons 
may occupy. The fourth and fifth step, the numerical 
solution of the Schrodinger equation, is a standard task 
that can be accomplished using solvers like the Lanczos 
or Davidson method in the time-independent case |39j . 
and propagators like the short iterative Lanczos [40] or 
the Crank-Nicolson [5T] method for the time-dependent 
version. For efficiency, they should be implemented using 
sparse matrix algebra. Note, however, that depending 
on the RAS constraints, the degree of sparsity may be 
smaller than in the FCI case. 

A true modification with respect to Full CI is made 
only in the second step, since RAS schemes require a 
more elaborate ordering and bookkeeping of the deter- 
minants. In order to construct the Hamiltonian matrix 
for the system ([2]), one needs to be able to calculate the 
matrix elements with Slater determinants \I),\J) 

Hjj{t) = {l\Hit)\j) = Y,hj,git){l\ala,\j) 

pq 

\a'^palasag\ j) . (14) 

pqrs 

In the last equality, we inserted the second quantization 
representation of the Hamiltonian, see e.g. Refs. ^25. .32] . 
which is given here as the sum of an explicitly time- 
dependent single-electron part (first sum containing ki- 
netic and potential energy and the external field) plus a 
two-particle part (the second with the Coulomb interac- 
tion). The operators Up (cij,) thereby annihilate (create) 
a particle in the orbital V'fc- The evaluation of the to- 
tal Hamiltonian is sketched in the following by means of 
the single-particle term. For a more detailed introduc- 
tion, we refer to Ref. [32] • First, one loops over the set 
of determinants | / ) and the orbital indices (p, q) [resp. 
(p, q, r, s)] and applies the excitation operators e.g. to 
the left, to obtain 

{i\d% = (15) 

with the phase factor ^ e {—1, 0, 1}. Note that one only 
gets a non-zero contribution if orbital p is occupied in 
I / ), and either orbital q is unoccupied ot p — q holds. 
Next, the term ^ hpq should be added to the Hamiltonian 
matrix element Hjj. One therefore needs an efficient 
scheme to retrieve the address add(| J)) of the deter- 
minant I J ) . The simple idea of a linear search in the 
determinant list is thereby inappropriate, as that would 
imply a quadratic scaling in the number of determinants. 



A more convenient method is to construct the addresses 
directly from the determinant by using graphical or com- 
binatorial techniques [211 li^Hii] . Following Ref. [IS] , in 
the Full CI case with M orbitals and N particles, one 
obtains 

N 

add(|/)) = 1 + J2zip,ip), (16) 
p=i 

with the function Z defined as 

m=A/-/+l L V /V / -I 

(17) 

for fc < and Z{N, I) = I- N. With this choice, the 
address function add(| / )) assigns a one-to-one mapping 
from the set of A'dct — (^) Slater determinants | / ) to 
the set of indices {1, • • • , A^dot}- 

As introduced above, RAS Hilbert spaces are com- 
posed as tensor products of a certain number of smaller 
FCI Hilbert spaces, and RAS determinants are given as 
tensor products of FCI determinants. It is thus not too 
difficult to evaluate the address in the FCI subspaces, 
and connect these values to obtain the address in the 
total RAS space. For a detailed derivation, see e.g. 

Refs. [lam]. 

III. APPLICATION TO PHOTOIONIZATION 

A. Spatial partitioning 

The RAS schemes applied in quantum chemistry are 
commonly considered in energy space, that is, the RAS 
constraints are chosen according to the orbital energies. 
Typically, a penalty is placed on the high-lying energy 
orbitals, since these are expected to contribute only little 
to the targeted groundstate or low excited states. This 
view is also indicated in Figs. [T] and [2| For an efficient 
treatment of photoionization processes, it is crucial to 
extend this concept also to the coordinate space, where 
the partitioning is, instead, performed according to the 
spatial regions. A basic example is shown in Fig. [3] 
There, the coordinate space is divided into a region in 
the vicinity of the nucleus (Si) and a region outside (S2). 
Over the region near the atomic core we construct the 
Hilbert space "Hi, which should appropriately describe 
the groundstate | ^"0 ) , while the Hilbert space H2 con- 
structed over the outside region is used to model the 
scattering states. In the same way as before, one can 
then restrict the allowed particle numbers; for instance, 
one could consider only single-ionization processes, and 
thus allow only for a single electron in B2 ■ Note that here 
one implicitly makes the assumption that the true singly- 
ionized states are reasonably well described by singly ex- 
cited determinants. In order to perform the spatial par- 
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FIG. 3. Spatial partitioning of the active space used in the 
RASCI treatment of photoionization processes. Bi marks the 
region in which the groundstate | ^'o ) is localized, B2 the con- 
tinuum. The example shows a four electron atom (e.g., beryl- 
lium) . 

titioning, a single-particle basis is required which is local- 
ized in the partitioned regions, as described in the next 
section. 

The introduced "division of space" -concept is tradi- 
tionally employed in the R-matrix method, which cur- 
rently presents the most successful approach to the pho- 
toionization problem of many-electron atoms. R-matrix 
calculations are usually performed in a time-independent 
formulation, e.g. by using Floquet theory [47], though 
recently also a time-dependent version has been pro- 
posed [30J. Shortly, the R-matrix method proceeds in 
the following way: After the division into an inner and 
an outer region, the wavefunction in the inner region is 
expanded in a set of angular-momentum eigenfunctions. 
The focus is mostly on single-ionization, i.e. only one 
electron may traverse to the outer region, where it is de- 
scribed, e.g., on a finite-difference grid |3S]; occasionally, 
also double ionization is considered |31| . The interac- 
tion of the outer electron with the inner electrons is de- 
scribed by using a multipole expansion, while the effect 
of the outer on the inner electrons is neglected. At the 
boundary, a set of equations has to be solved in order 
to connect the wavefunctions in the two regions, which 
are known as the R-matrix equations. Hence, the R- 
matrix method takes a different strategy than the RAS 
approach: while the R-matrix method starts from two 
separated spaces and connects them by a fitting proce- 
dure, the RAS scheme begins with a large connected Full 
CI space and retains only a selected subspace. 

In the formulation as stated before, the R-matrix 
method can be considered as a special case of the TD- 
RASCI method presented in this work (one could, how- 
ever, also think of more general expansions in R-matrix 
theory). In order to arrive there, one needs to apply the 
same spatial partitioning and allow only for single elec- 
trons in the outer region, while the wavefunction in the 
inner region needs to be expressed in angular-momentum 
eigenfunctions. Next, one has to approximate the elec- 
tronic Coulomb interaction accordingly, i.e. apply a 
multipole expansion to certain contributions resp. ne- 
glect other contributions. This would surely considerably 



speedup of the calculation. We note, however, that in its 
natural formulation the RAS scheme treats all the aris- 
ing interactions in an exact way, i.e. without any multi- 
pole expansion. Moreover, the RAS scheme offers several 
possibilities for extensions, which are simply achieved by 
changing the input parameters. For example, it may be 
easily extended to 

(i) the treatment of double ionization. This is easily 
accomplished by allowing for two electrons in the 
outer region. 

(ii) use more spatial subdivisions. For example, one 
could allow for double ionization only in a small 
region and for single ionization in a much larger 
region. 

(iii) apply a combined spatial and energy partitioning. 
For instance, the double ionization could be re- 
stricted to states with angular momentum I — 1, 
while in the other continuum states only single- 
ionization is allowed. This is useful, e.g., in the 
treatment of laser-assisted Auger decay. 

Note that for each of the mentioned scenarios, the 
equation of motion is simply given by the discretized 
TDSE and it is not necessary to derive specialized 
equations. The only task is to adjust the RAS restric- 
tions once in the beginning to the physical problem at 
hand. 

The R-matrix method had great success in the de- 
scription photoionization processes and provided several 
benchmark results for many-electron atoms. This makes 
us confident that the RAS scheme may be suitable as 
well, and that it hopefully can extend the range of appli- 
cations to scenarios which could not be treated efficiently 
so far. 



B. Spherical FEDVR basis 

As mentioned in the previous section, to apply the spa- 
tial partitioning the single-particle orbitals should be lo- 
calized in the spatial partitions. For the present treat- 
ment of photoionization, the basis should further be suf- 
ficiently flexible to describe the scattering states. These 
two requirements rule out several commonly used basis 
sets; for instance. Slater- and Gaussian-type functions, 
which are routinely used in quantum chemistry (and oc- 
casionally also in photoionization studies P5]), are not 
able to provide an adequate continuum. Other basis sets 
like Sturmian functions, though useful for photoioniza- 
tion [1], are not spatially localized. Here, we apply a 
partial wave expansion [25j . 

V'feim(r,ms) = r/™(0,</))a(ms), (18) 

r 

where denotes a spherical harmonic, and a a spin- 
eigenfunction {a € {a,/?}). In order to be localized. 
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the radial basisfunctions should be represented by finite- 
difference grids, B-splines or discrete variable representa- 
tions. Note that this choice implies that we can construct 
the partitions only with respect to the radial distance to 
the nucleus, but not to the angles. The use of orbitals 
of this kind further ensures that the Slater determinants 
are eigenfunctions of the orbital angular-momentum op- 
erator Lz and the spin-projection Sz- 

In this work, we use the finite-element discrete vari- 
able representation (FEDVR) basis for the radial basis- 
functions x(r) [50], which is illustrated in Fig. Q. The 
coordinate space is thereby divided into a chosen number 
of finite elements, and in each element the basisfunctions 
are given by (normalized) Legendre interpolating poly- 
nomials. 
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FIG. 4. Radial part of the single-particle basis ( 18 1 used 



which are constructed over a Gauss- Lobatto grid {rk} 
with integration weights {wk}- Additionally, at the 
boundary between two finite elements a bridge function 
is introduced to ensure the continuity of the wavefunc- 
tion. For a detailed construction of the FEDVR we refer 
to Refs. [50, 51 . For our purpose, the use of a FEDVR 
basis is of great convenience due to several reasons: first, 
the interpolating polynomials ( 19 ) are sufficiently flexi- 



ble to represent very general functions. Next, the matrix 
representation of local operators and the kinetic energy 
is sparse, which in turn leads to a sparse representation 
of the many-body Hamiltonian. And third, the finite- 
elements can be easily adjusted to model the chosen spa- 
tial division. 



C. Mixed basis 

In order to obtain a well prepared groundstate in the 
RAS calculations, one needs to perform a further ad- 
justment of the single-particle basis introduced before. 
Let us first shortly sketch the problem: Full Configura- 
tion Interaction Hilbert spaces are invariant under uni- 
tary rotations of the single particle basis. That is, the 
Hilbert spaces spanned by all iV-particle determinants 
constructed from a given single-particle basis {| V'fe )} of 
size Nil is exactly the same as if we would use any uni- 
tarily transformed basis | -(/'j ) = Ujk | V'fc ) • While be- 
ing mathematically equivalent, different choices in gen- 
eral differ in their numerical behaviour, and it is gen- 
erally advantageous to use basis sets which produce an 
as sparse representation of the Hamiltonian as possible, 
since this gives an efficient evaluation of matrix-vector 
products. In contrast, restricted active space calcula- 
tions are not invariant under unitary transformations of 
the single-particle basis. Most obviously, this can be seen 
for a single determinant, where the choice of the orbitals 
is of primary importance. As it is well known, the opti- 
mal set of orbitals for a single determinant (in the sense 



in this work, which consists of a combination of Hartree- 
Fock orbitals in the vicinity of the atom (left element), and 
FEDVR basisfunctions in the continuum (two elements on the 
right). In each element, the basisfunctions are constructed on 
a Gauss-Lobatto grid with 20 points. In the first element, the 
Is and 2s orbitals are emphasized. 



of giving the lowest energy) is obtained by solving the 
Hartree-Fock equations. This demonstrates a principle 
that is generally valid in RAS calculations: the more the 
RAS space deviates from the Full CI space, the better the 
single-particle basis needs to be adapted. Therefore, to 
be confident, RAS calculations are usually carried out 
in a well adapted basis, regardless of the accuracy of 
the many-body space. Typical choices are the Hartree- 
Fock basis [57] or the basis obtained from a Multiconfig- 
urational Hartree-Fock calculation. For both tasks, we 
employ our recently implemented Multiconfigurational 
time-dependent Hartree-Fock program [5S]. 

However, the usage of these adapted basis sets meets 
its limits in the treatment of photoionization. If the size 
of the basis Nj, becomes large, the creation and storage 
of the two-electron integrals [the gpqrs in Eq. ( 14 )] in 



the Hartree-Fock basis is a very tedious task, sinceboth 
require at least an effort of 0{N^). Thus, the advan- 
tages of the grid-like FEDVR treatment [which roughly 
scales with 0{N'^)] were completely lost. To circumvent 
this, we apply a mixed basis, in which we perform the 
transformation onto the Hartree-Fock orbitals only in a 
small region Bi around the nucleus, while the remaining 
part B2 is described by the grid-like basis (for an illus- 
tration see again Fig.js]). By this, we obtain both a good 
description of the groundstate, which is assumed to lie 
almost entirely in "Hi, and of the continuum, and this 
at a moderate effort. For the construction of the mixed 
basis, we first solve the (Multiconfigurational) Hartree- 
Fock equations in Bi using the spherical FEDVR basis, 
which provides optimized occupied orbitals. The virtual 
orbitals, however, are usually delocalized and do not re- 
semble excited atomic orbitals. Thus, to obtain appro- 
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priate pseudo-orbitals, we start from the ideal states and 
orthonormalize them to the bound states (for a more 
elaborate construction scheme, see [SS])- As a result, 
we obtain the unitary transformation matrix Ui, which 
transforms from the spherical FEDVR basis to the opti- 
mized orbitals. The matrix 

u - (-' ;) m 

is then used to transform the electron integrals in the 
total basis B, Eq. [7j The arising basis is illustrated 
in Fig. Q. A similar basis consisting of a mixture of 
Gaussian-type and FEDVR orbitals has been used in 
Ref. j6j to describe double ionization of Beryllium. 
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IV. RESULTS 

In the present section, we apply the TD-RASCI 
method to the atoms Helium and Beryllium, for which 
we calculate the photoionization cross sections from an 
explicitly time-dependent treatment. The present re- 
sults could alternatively be calculated using a variety 
of (mostly time- independent) methods like perturbation 
theory or Floquet theory, which are applicable to long 
pulse duration and yield sharp spectra. The following cal- 
culations can therefore be viewed as a proof-of-principle 
application, which constitutes a necessary first step for 
future application to explicitly time- dependent processes. 
Our main objective is thereby to check whether the rel- 
evant states exist and are located at correct energy po- 
sitions, so that they can participate in the simulations. 
Furthermore, we will not exploit the full computational 
power of the RAS scheme, but only use the two-fold di- 
vision of space illustrated in Fig. [Sj which is routinely 
applied in R-matrix calculations. Thereby, we divide the 
coordinate space in two regions, one in the vicinity of the 
atom (region 1) and one outside (region 2). The bound- 
ary is placed at i? = 20 bohr, which appears to be a 
common choice also used in Ref. [30]. In the outer re- 
gion, we allow for a single electron, i.e. we investigate 
only single ionization, while the groundstate in the inner 
region is described using different approximation levels. 
The basic approximation is TD-CIS, where the wavefunc- 
tion is described by one groundstate determinant plus all 
single-excitations. By including more determinants, we 
can then improve the description of the wavefunction in 
"Hi, and this way include relevant ionization channels. 

A. Helium 

The Helium atom is a thoroughly studied system for 
more than 40 years, and especially the last decade has 
seen an ever growing number of investigations, as, e.g., 
of pump- probe processes (53l IM]. Another topic of re- 
cent interest is the two-photon double ionization [T] [2] , 
which has also been investigated in our former work 



FIG. 5. (Color online) Single-ionization cross sections of He- 
lium, calculated for a laser field with intensity / = lO^^W/cm^ 
and a squared-sine envelope of duration T — 400 a. u. The 
black dots show experimental results of Samson [58], red 
squares the TD-FCI reference results and blue triangles the 
TD-RASCI results, which all agree almost perfectly. The re- 
sults of TD-CIS calculations are depicted by green triangles. 



using the Multiconfigurational time-dependent Hartree- 
Fock method 155]. Here, we use the TD-RASCI method 
to consider single photoionization for photon energies in 
the direct and sequential regime of two-photon double 
ionization. In Helium, the idea to partition the single- 
particle basis has been applied several times before, see, 
e.g., Ref. [55] for a spatial partitioning, or Ref. [2 for an 
optional partitioning in terms of angular momenta. Fur- 
thermore, R-matrix calculations have been performed on 
the single- and double ionization of Helium [ST] [56l [57] . 

To obtain appropriate Helium orbitals, we solved the 
Hartree-Fock equations in the region r < 20 bohr us- 
ing a radial basis consisting of 54 FEDVR functions for 
the Is orbital. For the virtual orbitals, we use single- 
hamiltonian eigcnfunctions up to /-symmetry, with the 
ideal ns orbitals orthonormalized to Is HF-orbital. For 
the TD-RASCI method, we then employed a number of 
10 s-orbitals and 3 p-orbitals to form the Full CI ground- 
state with vanishing spin- and angular-momentum pro- 
jection values. A single electron can be excited out of 
this space either into the remaining ideal orbitals inside 
r < 20 bohr or into the outside region, where a FEDVR 
basis is used with finite elements of length 4 bohr con- 
taining 10 basisfunctions each. The groundstate energy 
we obtain in this way is £' = —2.8894 Hartree (the ex- 
act non-relativistic energy is E — —2.9037 Ha). We have 
further applied orbitals which have been optimized with 
the MCTDHF method. Though they yield an improved 
groundstate energy of E = —2.9002 Ha, the photoioniza- 
tion spectra were less accurate. This is probably due to 
a poorer representation of excited states, which is caused 
by the fact that the orbitals are optimized for the ground- 
state. On the groundstate, we apply a laser pulse with 
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FIG. 6. (Color online) Single-ionization cross sections of He- 
lium for the same field parameters as in Fig. [5| B lack squares 
are the experimental results of Samson et aZ. p8j , and red ar- 
rows mark the resonance positions of ^P°-symmetry accord- 
ing to Scrinzi and Piraux 



FIG. 7. (Color online) Single-ionization cross section for 
Beryllium as a function of photon energy. Black squares de- 
note the experimental results of Wehlitz et al. [621 163| , the 
dashed curve the time-dependent calculations of Laulan and 
Bachau 5 . Both are compared to TD-CIS (green) and TD- 
RASCI (blue) calculations. 



a squared-sine envelope of duration T = 400 a.u. (~ 10 
fs) and intensity lO^^W/cm^. The ionization yields are 
calculated as the norm of the wavefunction in the region 
r > 20 bohr and are extracted at the end of the pulse. 
The yields are related to the single-ionization cross sec- 
tions through the formula given in Ref. [T]. 

Figure [5] depicts the cross section for photon energies 
in the direct regime of two-photon double ionization, i.e. 
w = 40 eV to w = 54 eV. Plotted are the results of TD- 
CIS, TD-RASCI, as well as experimental results of Sam- 
son et al. 58J. Further, we show time-dependent Full 
Configuration Interaction (TD-FCI) results taken from 
our former work j25j (calculated for a field duration of 
T = 100 a.u.). All curves agree qualitatively and show a 
smooth, monotonically decreasing behaviour. The TD- 
CIS curve, however, differs visibly from the experimental 
results, whereas the TD-FCI result, i.e. the exact solu- 
tion of the two-particle problem, agrees well with the ex- 
periments. The TD-RASCI results matches the FCI ref- 
erence result almost perfectly, which is remarkable as the 
calculations can be estimated to be about twenty times 
faster for the present parameters. One can thus conclude 
that, as expected, the largely reduced RAS Hilbert space 
suffices to describe the single-ionization, or, vice versa, 
that the excluded part of the Full CI space has a neg- 
ligible occupation. The lower accuracy of the TD-CIS 
result, on the other hand, demonstrates the importance 
of an accurate description of the electronic structure. 

In Fig. |6) we plot the cross section against photon 
energies in the sequential two-photon double-ionization 
regime, a; > 54 eV. At these photon energies, one ob- 
serves several Fano resonances in the cross-section, which 
correspond to ionization channels where doubly-excited 
states become occupied and subsequently decay via au- 
toionization occurs. The resonances arise due to the in- 



terference between these channels and the direct path- 
way, in which the electron is directly ionized into the 
continuum. The energies of the resonances with ^P°- 
symmetry are indicated by red arrows, and are taken 
from a work of Scrinzi and Piraux who applied the 
complex-scaling method [59j . The experimental results 
of Samson et al. i58j are denoted by black squares. They 
clearly resolve the resonance at 60.1 eV below the N = 2 
threshold, and the resonance at 69.9 eV below the N = 3 
threshold (for a more detailed classification we refer to 
|60j). The TD-CIS results yield a smoothly decreasing 
curve that shows no peaks. This is not surprising as the 
doubly-excited bound states are missing in the expansion 
of the wavefunction, so that only the linear photoion- 
ization process is possible. In contrast, the TD-RASCI 
wavefunction includes the autoionizing states and thus 
contains the relevant ionization channels by construction. 
The TD-RASCI curve agrees well with the experimental 
results, and shows peaks at the correct resonance po- 
sitions. Due to the restricted propagation time in our 
time-dependent method, the widths of the resonances are 
described less accurately, and much longer propagation 
times would be necessary to improve on that. For in- 
stance, the width of the 2s3p_ resonance requires a pulse 
of at least the same bandwidth, and hence a duration of 
at least 6 ps [61]. We stress, however, that the focus of 
the present method is not on calculating accurate cross 
sections, but rather to allow for a time-dependent treat- 
ment of the many-body photoionization problem. 
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FIG. 8. Single-photoionization cross section for Beryllium 
between the Be+(2p) and Be+(3s) thresholds. TD-CIS and 
TD-RASCI calculations are compared to the results of Kim et 
al. [65], which were obtained with the R-matrix method. 



B. Beryllium 

Beryllium is the simplest atom with two closed shells. 
The ionization energy of its 2s orbital is experimen- 
tally found to be 9.32 eV, and 123.35 eV for ioniza- 
tion from the Is orbital [M]. Therefore, when consid- 
ering photon-energies below the Is-threshold, the large 
energetic separation of the two shells allows for a fixa- 
tion of the Is core-orbitals and thus for the reduction 
to an effectively two-particle problem. In this spirit, 
several works have been performed focusing on single- 
and double-electron photoionization. Most of them are 
following the time-independent approach, using, e.g., 
the R-matrix method |65| , the relativistic random-phase 
approximation [M], hyperspherical calculations [H7], or 
the Multiconfigurational Hartree-Fock method [68]. The 
time-dependent treatments are usually performed with 
the close-coupling approach [31 [S] or a mixed-basis ap- 
proach [S] , both of which apply a reduction by freezing 
the Is orbitals (or, occasionally, the 2s orbitals). 

Here, for the TD-RASCI approach, we use a similar 
ansatz as before in Helium: after solving the Hartree- 
Fock equations and orthogonalizing the remaining ideal 
orbitals, we include 10 s-orbitals and 2 p-orbitals, which 
may be arbitrarily occupied by the 4 electrons. The only 
restriction is that we allow only for single-excitations 
from the Is orbitals, which is reasonable, as for the con- 
sidered photon energy range and intensity double excita- 
tions from the Is shell are negligible. By this, we obtain 
a RASCI groundstate energy oi E — —14.618 Ha, which 
has to be compared to the HF energy oi E = —14.573 
Ha and to the exact energy E — —14.667 Ha |69j. 

Figure [7] shows the cross-sections for photon energies 
in the range a; = 20 eV to 85 eV. The duration of the 
laser field is 400 a.u. and its intensity is 10^^ W/cm^. 



FIG. 9. Photoionization cross-section of Beryllium around the 
Is threshold. The red dashed line shows the results of Laulan 
and Bachau 5 obtained through a time-dependent method 
and fixation of the 2s orbitals, the green and blue curves our 
TD-RASCI and TD-CIS results. The black lines mark the 
resonance and ionization energy due to the NIST database. 



The experimental results (black squares) are taken from 
two works of Wehlitzei al. [551 [55]. Further, we show the 
theoretical cross section of Laulan and Bachau, which 
was calculated with a time-dependent method and fixed 
Is orbitals [5] as well. The TD-CIS results already show 
a good qualitative agreement. By improving the quality 
of the wavefunction, the TD-RASCI method is able to 
reproduce the experimental results almost perfectly. In 
Fig.jsj we consider the resonances in between the Be'''(2p) 
and Be''' (3s) threshold, which were obtained from a pulse 
duration of T = 800 a.u. The reference results are taken 
from an R-matrix study of Kim et al. |65j . One notices 
again, that TD-CIS by definition fails to model the peaks 
entirely, whereas the TD-RASCI method determines the 
gross structure correctly. Due to the limited bandwidth, 
our time-dependent approach does, however, not resolve 
the resonance region around 20 eV. To obtain more ac- 
curate results, much longer propagation times and also a 
more accurate RASCI expansion allowing for double ex- 
citations to d-orbitals would be necessary. Finally, Fig. [9] 
concentrates on photon energies in between w = 80 and 
150 eV, i.e. around the onset of the ionization of the Is 
orbital at 123.35 eV. Here, a pulse duration of T = 200 
a.u. Again, we compare our results to those of Laulan 
and Bachau, where this time the two Is electrons are 
considered as active, while the 2s shell is held frozen. In 
contrast. The TD-RASCI approximation (and TD-CIS 
as a special case of it) treats ionization from both shells 
on the same footing, and is therefore able to adequately 
describe the transition region. In TD-CIS, the ionization 
of the Is orbitals occurs at an energy of w = 128.72 eV, 
which is the ionization energy obtained in Hartree-Fock 
approximation through Koopmans theorem. The more 
accurate TD-RASCI wavefunction is able to shift the on- 
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set energy towards the exact position. The same holds 
for the Is 2s^ 2p -resonance, which, unhke the resonances 
considered so far, is ah'eady contained in the TD-CIS 
approximation, as the state is obtained from the ground- 
state by a single excitation from Is to 2p. Again, the 
TD-RASCI approximation corrects the resonance posi- 
tion, and almost achieves the exact result. Above the Is 
threshold, one encounters a large number of dense-lying 
resonances caused by the further excitation of the elec- 
trons from the 2s shell, which have been studied in de- 
tail by Voky et al. using the time-independent R-matrix 
method [70]. In contrast to the results of Laulan and 
Bachau, which show a smooth monotonically decreasing 
curve, the TD-RASCI method again reproduces the oc- 
curence and position of the resonances which get occu- 
pied in the shake-up process. 

V. CONCLUSION 

In this work, we introduced the time-dependent re- 
stricted active space Configuration Interaction (TD- 
RASCI) method for the ab-initio simulation of photoion- 
ization processes. Though well known in quantum chem- 
istry and based on a conceptually simple idea, to our 
knowledge the TD-RASCI method has in its full gener- 
ality not been applied to solutions of the time-dependent 
Schrodinger equation so far. Only specialized variants 
have been employed, like the single- and two-active elec- 
tron approximation or time-dependent Configuration In- 
teraction singles. In particular, the TD-RASCI method 
bears several similarities to the time-dependent R-matrix 
method, which can be considered the most successful ab- 
initio approach to many-body atoms at present. How- 
ever, though the underlying idea is nearly identical, the 
actual implementation differs considerably and is, in our 
opinion, conceptually much simpler for the TD-RASCI 
scheme. In fact, the only task is the selection of an ap- 
propriate active space, for which we presented a system- 
atic strategy. This accomplished, the derivation of the 
equations of motion proceeds as in the Full CI method. 



and is as easy as going from Eq. ([s]) to Eq. Q: all 
one has to do is to insert the restricted ansatz into the 
Schrodinger equation, and arrive at a matrix equation 
for the expansion cocfhcients. We further presented the 
specializations required to efhciently treat photoioniza- 
tion processes, namely the use of a spherical FEDVR 
single-particle basis and, in order to obtain an appropri- 
ate initial state, the transformation onto a mixed basis 
set. 

The TD-RASCI method has been applied to the cal- 
culation of single-ionization cross sections of Helium and 
Beryllium, and we particularly examined the quality of 
different RAS approximations. It was shown that time- 
dependent Configuration Interaction singles (TD-CIS) 
method is able to qualitatively describe the experimen- 
tally measured cross-sections in the linear photoioniza- 
tion regimes. However, by definition, it lacks the descrip- 
tion of resonances, simply because the relevant doubly- 
excited states are not included in the wavefunction. By 
taking them into account in the TD-RASCI method, we 
are able to appropriately resolve the resonances. Due 
to the time-dependent approach, the energy-resolution 
of the spectra is less accurate than in time-independent 
calculations. Yet, this is of no practical significance in 
time-dependent applications, where it only matters that 
the relevant states exist and are located at the correct 
energy positions. This requirement being fulfilled, we 
are confident that the TD-RASCI method is well suited 
for studying a large class of essentially time-dependent 
physical processes, including pump-probe scenarios or 
the Auger decay of many-electron systems. Our future 
work will therefore concentrate on an application of the 
TD-RASCI method to these problems. 
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